Anharmonic stabilization of the high-pressure simple cubic phase of calcium 
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The phonon spectrum of the high-pressure simple cubic phase of calcium, in the harmonic approx- 
imation, shows imaginary branches that make it mechanically unstable. In this letter, the phonon 
spectrum is recalculated using density-functional theory (DFT) ab initio methods fully including 
anharmonic effects up to fourth order at 50 GPa. Considering that perturbation theory cannot 
be employed with imaginary harmonic frequencies, a variational procedure based on the Gibbs- 
Bogoliubov inequality is used to estimate the renormalized phonon frequencies. The results show 
that strong quantum anharmonic effects make the imaginary phonons become positive even at zero 
temperature so that the simple cubic phase becomes mechanically stable, as experiments suggest. 
Moreover, our calculations find a superconducting Tc in agreement with experiments and predict an 
anomalous behavior of the specific heat. 



PACS numbers: 63.20.kg,63.20.dk,63.20.Ry,74.25.Kc 

The understanding of crystal lattice vibrations in 
terms of phonons provides an excellent paradigm to in- 
terpret and understand a wide range of physical phe- 
nomena [I . In most cases, the harmonic approximation 
describes accurately phonon frequencies and the associ- 
ated physical properties. However, there are examples 
where experimentally confirmed structures display imag- 
inary phonons in ah initio DFT calculations, indicating 
that in such cases anharmonicity cannot be neglected. 
The high-pressure simple cubic (sc) phase of calcium is 
an important example of the possible stabilizing role of 
anharmonicity. Indeed, while measurements confirm the 
presence and stability of this structure [2j-i5j , theoretical 
calculations based on the harmonic approximation find 
imaginary phonon branches throughout the whole Bril- 
louin zone (BZ) [6 9J. 

Under pressure, calcium exhibits a complex and inter- 
esting behavior. For instance, it becomes the element 
with the largest superconducting critical temperature 
(Tc), reaching 25 K at 161 GPa [T^. According to room 
temperature x-ray diffraction measurements [2l[3lfTl|[T2]. 
the ambient condition face-centered-cubic (fee) phase 
transforms to body-centered-cubic (bcc) at 19 GPa, to 
sc at 32 GPa, to P432i2 at 119 GPa, to Cmca at 143 
GPa and to Pnma at 158 GPa. Moreover, it has re- 
cently been reported that upon cooling the sc structure 
transforms into a very similar monoclinic Cmmm phase 
at 30 K and 45 GPa fS]. On the other hand, evolutionary 
DFT simulations within the generalized gradient approx- 
imation (GGA) at K [13] found that the experimental 
phases do not always coincide with the lowest enthalpy 
structures. This is quite dramatic in the stability range 
of the sc phase considering that the lAi/amd structures 
(from 33 to 71 GPa) and C2/c structures (from 71 to 89 



GPa) have considerably lower enthalpy than sc. Recent 
diffusion quantum Monte Carlo calculations (DMC) [9] 
have brought new insight to this problem, showing that 
the sc phase is energetically preferred over the lAi/amd 
phase when the exchange-correlation energy is treated 
correctly. Nevertheless, the question of dynamical stabil- 
ity remains and a proper quantum-mechanical treatment 
explicitly incorporating anharmonicity is still missing. 

The extreme anharmonicity in sc Ca requires a non- 
perturbative approach and suggests the application of 
the self-consistent harmonic approximation (SCHA) [TU 
Hg. The SCHA seeks the physically well-defined Gibbs- 
Bogoliubov bound and, in contrast to classical molecu- 
lar dynamics (MD) or statistical sampling methods [16], 
works at any temperature with no additional cost. How- 
ever, in order to apply this theory, the knowledge of all 
anharmonic coefficients is needed. Calculating them from 
first principles is usually complicated and highly time- 
demanding, thus, the SCHA has been normally applied 
making use of empirical potentials. Nevertheless, given 
the simplicity and high-symmetry of the sc structure, we 
have calculated ab initio all the necessary anharmonic co- 
efficients up to fourth order in displacement. The SCHA 
could then be applied to compute the temperature depen- 
dent renormalized phonon frequencies. The calculations 
have been performed at 50 GPa and, as it turns out, 
within this formalism the phonons of sc Ca are stable 
even at K at this pressure. Unless stated otherwise, we 
use atomic units throughout, i.e., h ~ 1. 

Within the Born-Oppenheimer approximation, the 
Hamiltonian describing the dynamics of the N ions in 
the crystal is given hy H — T '\-lJ , where T and U are, 
respectively, the kinetic and potential energy operators 
of the ions. The potential is expanded up to fourth order 
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a,sU = Uo + U2 + U3 + U4 with 

Un^^^ V ^i"^(qi)...?i""(q„)<i>"-""(qi,...,q„). 

{aq} 

(1) 

In Eq. ([T]), {a} represent Cartesian coordinates, 
$"i - ""(qi, . . . , q„) is the Fourier transform of the 71th 
derivative of the total energy with respect to the ionic 
displacements and (q) is the Fourier transform of the 
ionic displacement operator. In the harmonic approxi- 
mation, neglecting the third and fourth order terms of 
the potential, the Hamiltonian is diagonalized in terms 
of phonons. The term C/3 + U^^ can be treated within per- 
turbation theory to correct the phonon frequencies and 
account for their finite lifetime [TT. However, in sc Ca 
the energy has no lower bound due to the imaginary fre- 
quencies obtained in the harmonic approximation and, 
therefore, one needs to treat anharmonicity beyond per- 
turbation theory. In the SCHA, one adds and subtracts 
to the potential a trial harmonic term that yields well 
defined real phonon frequencies, C/j, and redefines the 
Hamiltonian as H = Hq + Hi , with Hq = T + t/" ^^^d 
Hi = (C72 - U^) + U3 + Ui- The exact free energy F 
satisfies the Gibbs-Bogoliubov inequality 

F<Fo + (i?i)o, (2) 

so that the minimum of the right-hand side of Eq. ([2]) 
becomes an excellent approximation of F . Fq and (-ffi)o 
are given as Fq = -^\nZ and {Hi)o = tr{Hie~'^"°)/Z , 
where /3 — l/fcsT and the partition function is Z = 

The adjustable parameters that can be used for the 
minimization are the trial phonon frequencies {Jlyq} that 
diagonalize Hq, v being a mode index. Differentiating 
Eq. ^ with respect to r2i/q, the equation for the trial 
frequencies that minimize the free energy can be obtained 
straightforwardly: 

^^q = ^Ici + ^^vciiu^- (3) 

A numerical solution of this equation leads to the renor- 
malized frequencies fJ^q at any temperature. In Eq. ([s]) 

X $"^"^"^'"*(q',-q',q,-q)[l + 2nB(l],,qO], (4) 

e^q is the phonon polarization vector, M the mass of Ca, 
{wi/q} the phonon frequencies diagonalizing i/2, imagi- 
nary at some q, and ub the usual bosonic occupation 
factor. As it can be seen, the third order anharmonic co- 
efficients do not contribute to F at this level of approx- 
imation. It should be noted that in the renormalization 
process the polarization vectors are kept fixed. This is 




FIG. 1. (Left panel) Harmonic phonon spectra and renormal- 
ized anharmonic phonon spectra at K and 300 K of sc Ca 
at 50 GPa. For the K anharmonic branches the value of 
the mode electron-phonon coupling, A^q, is proportional to 
the area of each filled circle. Filled squares depict the renor- 
malized frequencies obtained by Teweldeberham et al. [9] from 
classical MD at 300 K and 45 GPa. (Right panel) At zero tem- 
perature, the anharmonic results for the integrated electron- 
phonon coupling parameter A(tj), the Eliashberg function 
a'^F{u!) and the phonon density of states (PDOS). 



justified for the highly-symmetric sc phase, but in sys- 
tems with different atoms in the unit cell polarization 
vectors may be used to minimize Eq. (|2]). 

The computationally most expensive part of the 
method described above is the ab initio calcu- 
lation of the fourth order anharmonic coefficients 
{$aia2a3a4(q'^ -q',q, -q)}. These can be obtained tak- 
ing numerical second derivatives of dynamical matrices 
calculated in supercells (see, for example, Ref. [TS': the 
method presented there was slightly extended, given that 
two linearly independent real displacements must be used 
to generate the necessary supercells at q points not at the 
BZ edge). Such dynamical matrices were obtained using 
density functional perturbation theory (DFPT) [19] as 
implemented in Quantum ESPRESSO [20] within the 
PBE-GGA ^21j and making use of a 10 electron ultrasoft 
pseudopotential with 3s, 3p and 4s states in the valence. 
A 30 Ry cutoff was used for the plane-wave basis and 
a 16 X 16 X 16 Monkhorst-Pack k mesh for the BZ in- 
tegrations. Phonon frequencies and anharmonic coeffi- 
cients were calculated ona4x4x4q grid T2] and the 
renormalized phonon dispersion curves were obtained by 
Fourier interpolation. 

The strong anharmonicity in this system stabilizes all 
the imaginary phonon branches even at K, as can be 
seen in Fig. [l] This is an extraordinary effect consid- 
ering that, normally, anharmonic stabilization of unsta- 
ble modes occurs with increasing temperature 16J. MD 
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FIG. 2. Total energy per atom when the atoms are displaced 
along the transverse mode at M (a) and R (b). The dots 
depict the ab initio total energies, the solid line is the fit to 
the quartic potential, the dashed line the harmonic contribu- 
tion and the dash-dotted line ^huj. The $4,iyq values obtained 
from the fit to the quartic potential, frozen phonon (fp) re- 
sult, and from the differentiation of the dynamical matrices 
in supercells are shown. 

simulations have suggested that sc Ca may be stable at 
room temperature [9l [23] . On the other hand, MD cal- 
culations cannot [13] predict the stabilization below the 
Debye temperature {&d ^ 120 K according to our calcu- 
lations) since, as we demonstrate, sc Ca is stabilized by 
purely quantum anharmonic effects at K. In particular, 
our results give 26.8 cm~^ and 2.6 cm~^ for the trans- 
verse modes, unstable in the harmonic case, at the X 
and M points respectively. Although low-energy trans- 
verse modes show the largest renormalization, longitu- 
dinal modes also suffer a considerable anharmonic cor- 
rection. As expected, phonon frequencies increase when 
temperature is raised. Concretely, the temperature de- 
pendence is very strong for the transverse mode at M 
and at R. At 300 K, above 9d, our results are in close 
agreement with the values obtained from MD at the zone 
boundary by Teweldeberham et al. [5]. Note that their 
calculation was performed at 45 GPa and ours at 50 GPa. 

In sc Ca, as in many other cases [13 HH Hi], it is 
crucial to account for scattering between phonons with 
different momenta. Indeed, if we make the assump- 
tion that nondiagonal coefficients are equal to the di- 
agonal ones in Eq. Q, $"i"2°3a4 (q' _q'^ q _q) _ 

<I>"i"2"3"*(q, — q, q, — q), the error caused in is quite 
dramatic and the temperature dependence becomes un- 
realistic. For example, the mode at R reaches a frequency 
of 235.5 cm-i at K and aheady 372.2 cm"! at 100 K, 
2.6 times larger than our predicted value. At R such a 
difference is a consequence of the huge value of the diag- 
onal coefficients in comparison to the nondiagonal ones. 
As shown in Fig. [2] the importance of the diagonal anhar- 
monic coefficients can be calculated from frozen phonon 
calculations. For a mode with momentum q at the edge 
of the BZ, when the atoms are displaced from their equi- 



librium position R as rya cos (qR)ei,q, with a the lattice 
parameter and rj a small dimensionless number, the po- 

2 4 

tential is given as E/N{r]) = \o?MijjI^ + ^$4,i/q, with 

{a} 

(5) 

A fit to this potential yields the frozen phonon $4,iyq coef- 
ficient. As can be seen in Fig. [2] our values obtained dif- 
ferentiating dynamical matrices in supercells are in good 
agreement with frozen phonon estimates. 

Our method yields the whole renormalized phonon 
spectrum at K and, thus, we can estimate the super- 
conducting transition temperature in sc Ca. The usual 
electron-phonon vertex is not modified by anharmonicity 
since the matrix elements of the gradient of the crys- 
tal potential are independent of the phonon frequen- 
cies [26|- Therefore, the electron-phonon coupling con- 
stant A can be calculated straightforwardly using the 
electron-phonon matrix elements and the renormalized 
frequencies fi^q at K. The convergence of the electron- 
phonon matrix elements required a denser 80 x 80 x 80 
k grid. Integrating the Eliashberg function, a'^F{ui), we 
obtain A = 0.74 and uJiog = 53 K, leading to Tc ~ 2.1 
K, an estimate obtained from the Allen-Dynes modified 
McMillan equation [57] (we have used /i* = 0.1 for the 
Coulomb pseudopotential). As can be seen in Fig. [ij the 
greatest contributions to A come from the soft transverse 
modes which are unstable in the harmonic approxima- 
tion. Indeed, the integrated electron-phonon coupling 
parameter, reaches the value of 0.54 at 50 cm~^ 

(that is, 73 % of the total value of A), so that if it were 
not for these soft modes sc Ca would superconduct only 
below 0.1 /xK. The value calculated for Tc at 50 GPa 
is in close agreement with the experimental 1.2 K value 
obtained by Okada et al. [55] and with the 1.7 K value ob- 
tained extrapolating linearly the Tc values measured for 
sc Ca at higher pressure in more recent experiments [10) . 
Finally, despite the strong anharmonicity, the isotope co- 
efficient, a — " ^}"^/ , is predicted to be 0.45, close to the 
0.5 value of a BCS superconductor. 

The temperature dependent renormalized frequencies 
{^^i/q} allow us to estimate the anharmonic free energy 
directly from Eq. ^ and the constant volume specific 

heat as Cy = —T ^fjS^ . As shown in Fig. jsj the high- 
temperature limit of the specific heat per atom is reduced 
by 17 % from the classical Zks value given by the equipar- 
tition theorem. Such a reduction from the Dulong-Petit 
law is a sign of strong anharmonicity [551 132] and has 
been observed in different systems [3T1 [35]. Moreover, 
the low-temperature behavior of Cv is strongly modified 
from the relation of harmonic crystals. The anomalies 
of the specific heat are mainly driven by the temperature 
dependence of the phonon frequencies in F^. Indeed, as 
depicted in Fig. [3] when the specific heat is calculated 
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FIG. 3. (a) Constant volume specific heat per atom for tlie 
sc Ca anfiarmonic crystal (line with grey circles). For com- 
parison, the specific heat calculated with the usual harmonic 
formula including the renormalized Jl^q frequencies at K is 
shown (solid line). The Sfcs line is depicted too (dashed line), 
(b) Low temperature specific heat in logarithmic scale. 



from Fq assuming that the K renormalized phonons 
are not modified under temperature, the Dulong-Petit 
law and the low-temperature behavior are recovered 
as expected. 

In summary, within the SCHA, using a variational pro- 
cedure based on the Gibbs-Bogoliubov inequality we have 
shown that the high-pressure sc phase of Ca is stabi- 
lized even at K by strong quantum anharmonic effects. 
This procedure has been used calculating fully ab ini- 
tio the anharmonic coefficients up to fourth order in the 
whole BZ and may be applied as well in many cases where 
the phonons are imaginary or anharmonicity needs to be 
treated beyond standard perturbation theory. Although 
below 30 K the sc phase may transform to a rather similar 
monoclinic Cmmm phase [S], which is mechanically un- 
stable in the harmonic approximation as well and shows 
very similar harmonic phonons [9], we have calculated 
the superconducting of sc Ca finding a good agree- 
ment with experiment. Moreover, the huge anharmonic- 
ity in this system makes the specific heat very anomalous 
with a large reduction from the high-temperature 3k b 
limit. An experimental confirmation of this last feature 
would indirectly show the strong anharmonic behavior 
predicted. 
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